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7 .STOCHASTIC  INVERSE  PROBLEM  IN  THE JIADIATION  OF  NOISE^ 

P l/cHOW*»«B»  L^/maESTRF.LI.C||I 

Abstract.  A general  procedure  for  regularizing  the  stochastic  inverse  source  problem  is  presented  The 
present  approach,  based  on  the  minimization  of  a structural  functional,  allows  the  incorporation  of  a priori 
knowledge  about  the  structure  of  the  undetermined  source.  As  special  cases,  it  yields  the  methods  of 
Lagrangian  multipliers.  Tihonov's  regularization  and  the  pseudoinverses.  Possible  generalization  and 
connection  with  Franklin's  well-posed  extensions  are  discussed  briefly 

1.  Introduction.  This  paper  is  concerned  with  an  inverse  problem  in  the  radiation 
irom  a randomly  fluctuating  source.  The  problem  was  suggested  by  a study  on  the 
determination  of  the  unknown  source  distribution  in  a jet  from  the  far-field  acoustic 
radiation  data.  According  to  Lighthill's  model  ( 1 j,  the  jet  noise  is  generated  by  the 
fluctuating  part  of  the  fluid  velocities  in  the  jet.  In  this  case,  the  pressure  field  is 
governed  by  an  inhomogeneous  wave  equation,  where  the  source  term  is  given  by  the 
second  order  divergence  of  the  Reynolds  stress  tensor.  A great  majority  of  works  [2] 
in  jet  noise  is  concerned  with  finding  the  pressure  correlation  function  by  specifying 
the  source  covariance  function,  e g.  a quadrupole  of  certain  form.  In  reality,  the 
source  is  neither  known  a priori,  nor  directly  measurable.  It  is  the  acoustic  pressure 
fluctuation  at  a given  point  that  can  be  measured  by  a microphone.  Question  naturally 
arises  as  to  the  possibility  of  using  a set  of  acoustic  radiation  data  to  infer  the  charac- 
teristics of  the  source  distribution.  This  gives  rise  to  a stochastic  inverse  problem. 

The  common  mathematical  difficulty  in  dealing  with  inverse  problems  in  physical 
sciences  lies  in  the  fact  that  they  are  generally  ill-posed  in  the  Hadamard  sense  [3|. 
Recall  that  a mathematical  problem  is  well-posed  if  the  problem  has  a solution,  the 
solution  is  unique  and  it  depends  continuously  on  data.  The  inverse  source  problem  is 
clearly  ill-posed,  because  a set  of  radiation  data,  say,  on  a surface  enclosing  the  source 
region  is  insufficient  to  determine  the  source  distribution.  However  the  partial  and 
inaccurate  data  should  be  of  value  in  the  reconstruction  of  source. 

Conceivably  there  are  many  possible  ways  to  utilize  the  radiation  data  to  yield 
some  information  about  the  source.  A procedure  that,  by  introducing  a priori 
assumptions  or  otherwise,  renders  a unique,  stable  solution  to  an  ill-posed  problem  is 
called  a regularization  or  a well-posed  extension.  Most  of  regularization  procedures 
clusters  around  the  well-known  method  of  least  squares,  notably  Tihonov’s 
regularization  |4]  the  well-posed  extension  by  Franklin  [5]  and  others  [6],  [7],  Though 
the  approaches  are  different,  these  methods  for  regularizing  deterministic  problems 
are  not  totally  unrelated.  The  present  problem  is  stochastic  in  nature  and,  hence,  there 
is  need  for  a generalization  of  the  current  procedures.  However  a full  generalization 
will  not  be  done  here.  For  trace  ability,  our  analysis  will  be  confined  to  the  estimation 
af  the  source  covariance,  known  as  the  correlation  theory  in  engineering  and  science. 
Then,  within  this  theory,  we  are  able  to  adapt  a generalized,  deterministic  method  to 
ihe  stochastic  problems.  The  regularization  procedures  of  Tihonov  and  others  are 
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the  physical  problem  should  reduce  the  degree  of  uncertainty,  we  shall  propose  a 
variational  regularization  procedure  which  allows  the  incorporation  of  a priori  knowl- 
edge about  the  structure  of  the  undetermined  source.  Our  procedure  resembles  the 
common  practice  of  introducing  an  extremal  principle  in  continuum  physics  to 
characterize  the  physically  acceptable  solution  among  admissable  ones.  To  this  end,  a 
structural  functional  and  a residual  functional  depending  on  a parametric  function  A 
will  be  introduced.  The  sum  of  these  two  will  be  termed  the  Lagrangian  functional. 
Then  our  regularization  procedure  consists  of  minimizing  the  Lagrangian  functional, 
with  or  without  constraints.  It  is  believed  that  the  proposed  method  is  more  flexible 
and  better  adapted  to  the  physical  problem  under  consideration. 

To  be  specific,  this  paper  will  be  concerned  with  a stochastic  inverse  radiation 
problem  in  a uniform  medium.  The  formulation  and  a preliminary  analysis  of  the 
inverse  problem  are  given  in  § 2.  As  an  example,  a simple  model  consisting  of  an  array 
of  point  sources  is  presented  and  analyzed  in  § 3.  Here  the  entropy  functional  is 
chosen  to  be  the  structural  functional  in  determining  the  source  distribution.  In  the 
following  section,  a general  theory  for  the  stochastic  inverse  problem  is  introduced. 
For  computational  feasibility,  our  emphasis  will  be  placed  on  quadratic  structural 
functionals.  It  will  be  shown  that,  by  specializing  the  structural  and  residual  functionals, 
the  general  procedure  yields  the  methods  of  Lagrangian  multiplier  [8),  Tihonov’s 
regularization  and  the  method  related  to  the  generalized  or  pseudoinverses  [9],  This  is 
done  in  § 5.  In  the  Appendix,  a possible  generalization  to  other  radiation  problems  and 
a modified  version  of  Franklin’s  stochastic  extension  are  briefly  discussed. 

Inverse  problems  arise  from  all  areas  of  engineering  and  physical  sciences.  The 
well-known  examples  are  the  determination  of  historic  climate  by  solving  the  heat 
equation  backward,  finding  the  mass  distribution  in  the  earth  by  measuring  the 
gravitational  potential,  the  inverse  scattering  problems  in  the  classical  and  quantum 
mechanics,  etc.  References  to  these  problems  can  be  found  in  the  books  [10]-[I2].  An 
informal  and  interesting  introduction  to  the  subject  of  inverse  problems  was  given  by 
Keller  [13].  The  nature  of  nonuniqueness  of  solutions  poses  a much  more  serious 
question  [14],  which  will  be  our  major  concern.  Therefore  a regularized  problem,  by 
convention  in  this  paper,  will  often  refer  to  a problem  with  a unique  solution.  The 
continuous  dependence  question  can  be  settled  separately  by  an  additional  smoothing 
procedure,  if  necessary. 

2.  Formulation  and  preliminaries.  Consider  the  acoustic  radiation  problem  in  a 
uniform  medium  due  to  a randomly  fluctuating  source.  Let  the  source  distribution 
q(x,  t,  w)  be  a random  function  defined  in  a region  D in  R'\  time  t SO  and  w 6 Q being 
the  underlying  sample  space.  Then  the  fluctuating  pressure  field  p(x,  I,  at)  satisfies  the 
inhomogeneous  wave  equation 

(2.1)  = -f€R’  />0’ 

where  the  speed  of  propogation  c is  constant  and  q = 0 if  x £ D.  In  a theory  of  jet 
noise,  q = (tf/bx,  ■ Hx k)T,k,  where  Tlk  is  known  as  the  Lighthill’s  turbulent  stress  tensor 
and  the  summation  convention  is  in  effect.  Given  q,  the  pressure  field  p in  (2.1)  is 
obtained  in  terms  of  the  retarded  potential  1 1 5 1 


, i f 


(r/c),u>)  - 
"4. 


(2.2) 

where  r - |jf  -f|. 
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For  a smooth  random  function  f(x,  t,  to),  its  ensemble  average  (or  mathematical 
expectation)  over  the  sample  space  fl  is  denoted  by  (/).  Assuming  the  source  q being 
time-stationary,  the  mean  fields  (q)  and,  hence,  ( p ) are  independent  of  time.  They  may 
be  taken  to  be  zeros  without  affecting  the  acoustic  signal  fluctuations.  Then  the 
covariance  of  p and  q become  simply 


Hx.  9,  r)  = (p(x,  t)p(y,  / + r)>, 
(2-4)  0(x,y,T)  = (q(x,t)q(y,t  + T)). 


In  view  of  (2.2),  the  following  integral  equation  can  be  derived 

M Hi. ,)- (±)2 1 _„<T<„ 


where  D2  = D x D,  r = |*  - £ |,  j = |y  - yj|. 

Now  let  Sp  be  a spherical  surface  enclosing  the  region  D with  radius  p.  Suppose 
that  the  radiation  data  P is  given  for  all  x,yeS„  and  |r|  < oo  as  p -*  oo,  (in  reality,  only  a 
finite,  discrete  set  of  data  is  obtainable).  One  is  required  to  determine  the  covariance 
function  <5  in  D based  on  an  incomplete  and  inaccurate  set  of  data  on  S = Sp  as  p -» oo. 
As  it  stands,  the  problem  is  ill-posed.  Aside  from  the  question  of  continuous  depen- 
dence, there  may  be  too  many  solutions. 

To  simplify  the  relation  (2.5)  in  the  far  field,  we  assume  that  the  spectral  density 
functions  P,  Q exist.  They  are  defined  as  the  Fourier  transforms  of  P,  0 in  r with  the 
parameter  (kc): 


(2  6)  f „ 

P(x,  y,  Ar ) = J P(x,y,T)e-'k‘'dT, 

(2  7)  0(x,  y,  k)=  | <5(f,  y,  r)e~'k‘' dr. 

Then  we  have  P*  (x,  y,  k ) = P(y,  x,k)  = P(x,  y,-k ) and  P(x,  x,k)=  P(x,  x,  -k ) g 0, 
where  means  the  complex  conjugate.  Similar  relations  hold  for  Q. 

Noting  (2.6)  and  (2.7),  with  a Fourier  transform  of  (2.5),  we  obtain 


(2.8)  P(i,y,k)=^  G(r,k)G*(s,k)Q(iTj,k)didv,  ^ 0 Q f 


where 

7 

?|f? 

(2-9)  G(r,k) 

e‘k'  l 

4 nr 

J 

l • 1 

is  the  free-space  Green’s  function. 

If  |f  | = |y  | = p,  as  p oo,  we  get  Vs^->  \ S 

(2. 10)  r = |jc  • £ 

s = \y-n\~P-0  v 

Here 

(2.1  I)  a=~. 

P 

p 

are  unit  vectors  along  x,  y,  respectively. 
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Upon  using  (2.10)  and  (2.11)  in  (2.8),  we  obtain 

P(x,  y,  k)~  G(p,  k )G*(p,  k)  f 0(0  v ,k)e  ‘kti  * ' n didfj , 

Jd2 

(2.12) 

x,  y € p - oo. 

Let  R(a,0,k)  denote  the  radiation  amplitude,  which  is  defined  as  the  limit  of  the 
normalized  covariance. 


(2.13) 


R(a,  0,  k)~  lim 


P(x.y.k) 


i.i*sG(p,k)G*(p.kj 

p-ac 


By  means  of  (2.13),  the  asymptotic  relation  (2. 12)  yields  an  integral  equation  for  R: 
(2.14)  R(a,0,k)  = f 0(i,  v,  k)e~'k<a  i') d0dr),  a,0eU, 

•I  D2 

where  U is  the  unit  sphere  centered  at  the  origin.  Let  R = Ri  + iR2  and  Q = Oi  + i02 
so  that  R , and  <?,  are  real,  / = 1,  2.  By  taking  the  real  and  imaginary  parts  of  (2.14), 
one  obtains  the  following  pair  of  coupled  integral  equations: 


(2.15) 


(2.16) 


R{(d,j3,k)=  k) cos  k(d  f-0  rj)didrj 

Jd 2 


-I 


02(lv,k) 


sin  k(a  ■ (-0  fj ) d\ dr\. 


R2(d , 0,  k)=  Oi(l  h.  *)sin  k(d  0-0  ■ rj)d0dr\ 

J D2 

+ [ 02(i<  V’  k)cos  k(a  ■ 0-0  rj)d£ dfj. 

Jd2 


We  note  that  the  integral  operator  in  (2.14)  takes,  say,  an  integrable  function  on 
D ‘ of  six  dimensions  into  an  integrable  (in  fact,  continuously  differentiable)  function 
of  the  product  sphere  t/2==t/xt/  of  a four-dimensional  manifold  Therefore  the 
integral  equation  (2.14),  or,  equivalently,  the  equations  (2.15)  and  (2.16),  merely 
specifies  a set  of  admissible  solutions,  if  it  exists,  instead  of  determining  a unique  one. 
The  question  of  how  to  provide  a logical  choice  of  sensible  solutions  will  be  our  main 
concern. 


3.  Analysis  of  a discrete  model.  Before  presenting  a general  theory,  it  is  instruc- 
tive to  work  out  a simple  model  in  detail.  Let  us  consider  the  discrete  model  of  N 
point-sources  located  at  x,  = (a,,  0, 0),  /'  = l,  2,  ■ ■ • , N,  on  the  x-axis,  where  x is  the  first 
coordinate  of  x.  Assuming  that  the  point  sources  are  isotrophic,  the  source  function  q 
takes  the  form 


(3.1) 


N 


q(x,l,to)=  X 

n I 


q,(l,to)fi(x  -x,). 


where  fi(x)  denotes  the  Dirac  delta  function  and  q,(t,  w)  arc  independent,  centered, 
second  order  (weakly)  stationary  random  processes.  Prom  (3.1 ),  we  form  the  covari- 
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ance  0 defined  by  (2.4) 

(3-2)  <5(x,y,r)  = £ (5,(r)«(x -x,)6(y 

i - 1 


with 

(3-3)  0,(r)  = (q,(t  )q,(t  + r)), 

or,  by  taking  the  Fourier  transform  (2.8), 


(3.4) 


<?(i.  M)=  I Ot(k)S(x  - x,)8(y  - x,). 

i-i 


Accession  i cr 

MIS  G^Ai'1659- 
K>C  TAB 
Unannounced 
Justification 

| alZL 

LDlst  -ibut  i / 

I Ave  i 7 " T-. ; . 1 j f y. 

;■  ■ ' ' " i/or 


A substitution  of  (3.4)  into  (2.16)  and  (2.17)  yields 
(3.5) 


N N 

/?i(y.  k)=  £ (?,(*)  cos  kya„  R2(y,k)  = £ (?,(<:)  sin  kyah 

/-•  /-i 


Here  y = cos  0O  - cos  and  0a,  are  the  angles  between  a,  0 and  the  x-axis, 
respectively,  and  a,  0,  x,  are  coplanar. 

Suppose  that  on  the  measurement  surface  S,  a finite  set  of  radiation  data 
{/?(»,  k),j  = 1,  2,  • • • , M)  is  taken.  Put 


(3.6) 


rm(k)=  Rx(ym,k),  rM+m(k)=  R2(ym,  k),  m = 1,  2,  • • • , M. 


em.n(k)  = coskyman 


eM.m.„(k)  = sin  ky„a„ 

n = 1,  2,  • • ■ , Af. 


m = 1,2,  •••,*/, 


Then  (3.5)  can  be  written  as 


(3.7) 


I em.n(k)Qm(k)  = r„(k),  m = 1, 2,  ••• , 2M. 


Let  E = (em„)  be  the  coefficient  matrix  and  f?  be  the  augmented  matrix  by  adjoining 
the  vector  (rm)  to  the  last  column  of  E.  It  is  well-known  [8]  that  (3.7)  has  a solution  if 
and  only  if  the  rank  r(E)  of  E is  equal  to  that  of  /?.  The  solution  is  unique  if  r(E)=  N, 
which  is  not  true  in  general.  For  definiteness,  we  assume  that  r(E)<N,  so  that  the 
linear  system  (3.7)  is  underdetermined.  To  single  out  a “reasonable”  solution  of  (3.7), 
we  may  apply  one  of  many  versions  of  the  least-square  method,  mentioned  in  the 
previous  section.  As  an  alternative,  we  introduce  a selection  principle  based  on  the 
entropy  function  9[0  1 of  noise  defined  by 


(3.8) 


nO]  = - I \ Qn(k)\nQ„(k)dk, 
n-i  J 


and  seek  among  those  solutions  of  (3.7)  the  one  that  maximizes  9 (or  minimizes  -9"). 
According  to  this  procedure,  we  regard  the  most  favorable  solution  of  (3.7)  as  being  the 
set  of  point  sources  that  are  distributed  in  the  most  chaotic  manner.  Here,  of  course,  we 
think  of  Ov’s  as  the  probability  density  functions  in  the  information  theory  (16).  Now 
the  regularized  problem  can  be  solved  by  the  method  of  Lagrangian  multipliers  [8],  To 
this  end,  let  \m(k),  m = 1, 2,  • • • , 2 M,  be  functions  defined  for  |j(|<ao,  and  let  us 
introduce  a residual  functional  #k\Q)  of  the  form 


*a|0|=  1 [*».(*)[  £ em.n(*)a(*)-Ym(*)]rf*. 
(n-i  J Ln  I J 


(3.9) 
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In  this  case,  \m  are  just  the  Lagrangian  multipliers,  but,  later  on,  they  may  be 
unspecified  parameters.  Let  if*  denote  the  Langrangian  (functional)  defined  as 

(3.10)  if*[0]  = ^[C>j+sr*[0]. 

For  a maximal  ;T,  its  first  variation  5if*  with  respect  to  Q must  vanish 

Nr  2 M N r 

(3.11)  S^*[<?]  = - I [1+ln  <?„]*?„  V £ kmem.„6Q„  dk  = 0, 

n-l  J m = I n * 1 * 


which  gives  the  Euler  (or  Lagrangian)  equation 

2M 

(3.12)  In  (?„(£)  = I em„(k)\m(k)~  1, 

m - 1 

or 

(3.13)  <?„(*)  = exp  | 2I  e„,„(k)Am(k)-l},  « = 1, 2,  • • • , N. 


The  above  formula  determines  the  source  correlations  O.'s  provided  that  A„’s  are 
known.  To  find  Am’s,  we  substitute  (3.13)  into  (3.7)  to  get 

(3.14)  I em.„(*)exp  ! I e,.„(k)Ai(k)- 1 1 = r„,  m = 1,  2,  • • • , 2M. 

* - 1 1 

The  above  system  of  nonlinear  equations  constitutes  2 M smooth  surfaces  in  A -space. 
If  these  surfaces  have  a unique  point  of  intersection  A = (A  t , • • • , A2*f),  then  the 
corrsponding  source  distribution  (5„  given  by  (3.13)  is  the  desired  solution.  As  a trivial 
example,  if  the  auto-correlation  r(k),  or  the  mean-square  intensity  of  radiation,  at  a 
given  point  on  S is  the  only  measurement  taken,  then  the  svste  ms  (3.7)  and  (3.13) 
reduce  to 

(3.15)  I 0 ,(*)=<■(*), 

n = I 

(3.16)  Q„(*)  = exp{A(*)-l},  n = 1,  2,  • • • ,N. 

In  this  case,  we  infer  that  the  point-sources  are  identically  distributed  with  Q„(k)  = 
r(k)/N,  n = \,2,-  -,N. 

In  view  of  (3.14),  a regularization  by  the  entropy  function  3~\0\  leads  to  a 
nonlinear  system  for  Am’s.  From  the  computational  viewpoint,  a linear  system  is  much 
preferred.  Therefore,  for  the  general  problem,  the  regularization  by  a quadratic- 
functional  & will  be  emphasized. 

4.  A general  theory  of  stochastic  inverse  source  problem.  In  general  we  must 
consider  a source  distribution  whose  covariance  <3(x,  y,  r)  is  a generalized  function 
Without  complicating  the  matter  by  invoking  the  theory  of  generalized  functions  [17], 
0 will  be  taken  to  be  a square-integrable  function  over  D’  xR,  This  amounts  to 
restricting  our  attention  to  almost  continuous  source  distributions.  Howes er.  since 
each  generalized  function  may  be  regarded  as  the  weak  limit  of  a sequence  of  ordinary 
functions,  the  result  in  this  section  will  still  be  valid  for  singular  source  distributions,  if 
interpreted  in  the  weak  sense. 

For  convenience,  set  y - («,  0).  i =■  (f,  f),  ( = (f.  fj)  and 

K(y,  z,  k)- exp-{ik(a  ■ (- 0 rj)).  KO(y,  k)^[  K(y,  (.  k )(?(<,'.  k ) 


(4.1) 
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Then  the  integral  equation  (2. 16)  can  be  written  as 

(4.2)  KQ(y,  k)  = R(y,  k),  yet/2. 

For  each  k,  the  integral  operator  K takes  a complex-valued  function  (?e  L2(D2)  into 
a complex-valued  function  R e L2(U2). 

For  the  purpose  of  regularization,  let  us  introduce  a structural  functional  ?T  of  the 

form 

(4.3)  J[Q,  O*)  = J J F[Q(z,  k ),  (?*(£,  k );  VQ(£,  k ),  VO*(z,  *);£,*]  dz  dk. 

Here  the  real-valued  function  F is  assumed  to  be  continuous  in  all  of  its  arguments, 
continuously  differentiable  with  respect  to  O,  Q*  and  their  gradients  VO,  VO*  in  z. 
The  functional  J plays  the  role  of  a constitutive  relation  in  continuum  physics  in  the 
form  of  an  extremal  principle.  The  functional  should  be  constructed,  insofar  as 
feasible,  by  taking  into  account  all  of  the  a priori  information  about  the  noise 
structure.  If  none  is  available,  a suitable  J can  be  chosen  for  the  purpose  of  smoothing. 
To  define  5",  let  Hi  be  the  subspace  of  H = L2(D2  x ft)  consisting  of  complex  functions 
O in  H,  whose  first  derivatives  dQ/dZj,  / = 1, 2,  • • • , 6,  are  also  in  H.  In  view  of  (4.3) 
and  the  ensuing  assumptions,  the  structural  functional  & is  well-defined  in  Hi  and  its 
first  variation  is  given  by  (18) 

8J[Q,  O*)  = | J {(F0  - V Fvo)«C?  + (F0.  - V • Fvo*)«0*>  dz  dk 

(4.4) 

+ f 

JaoJ 


{Fvq80  +F?o’0*)  ' n dz  dk< 


where  Fq  - dF/fiQ,  Fvo  stands  for  the  gradient  of  F with  respect  to  the  complex 
vector  VO,  and  dD 2 is  the  boundary  of  D 2 with  the  outward  normal  n.  Since  F is  real, 
we  have  Fq  — Fq-  and  F?0  =Fvo‘  Now  we  let  E(u)  be  a smooth  function  of  a 
complex  variable  u such  that  |E(«)|  is  convex  and  E(0)  = 0,  e.g.  E(u)  = 
u,  uu* , • , etc.  For  every  complex  A(y ,k)e  L2(U2  x R),  the  residual  functional  is 
defined  as 


*UQ)=  {A(y,*)F[KC>(y,*)-K(y,*)) 

(4.5)  VJ 

+ A*(y,  k)E*[KO(y,  k)-R(y,  k)])dydk, 

which,  by  construction,  is  real-valued. 

As  in  the  discrete  model,  we  introduce  the  Lagrangian  functional 

(4.6)  2>AQ)=‘r\0,0*]  + %*\0], 

but  the  parametric  function  A need  not  be  interpreted  as  a Lagrangian  multiplier.  In 
the  regularization  process,  wc  propose  the  principle  of  extremal  Lagrangian  with  or 
without  side  conditions.  To  proceed,  the  existence  of  extrema  will  always  be  assumed. 
For  a stationary  yA,  we  get,  by  invoking  (4.4M4.6) 

&y\|0|-S.T(a0*|  + *#A|0|  0, 


(4.7) 
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or  < 

j J {(Fo  - V • Fro)5Q  + ( Fo - - V F*0-)5(?*}  dz  dk  j 

(4.8)  + [ f [ {AKF(KQ-R)5Q+A*tf*£*(KQ-F)5Q*)di  dy  dk 

JD2  Ju2  J J 

+ f [ (Fvq  ■ hSO  — Fvq*  ' hSO*)  dz  dk  — 0, 

hD2  J j 

in  which  £(u)=  dE(u)/dul  and  M|  = Re{w}.  Since  SO  and  SO*  are  two  arbitrary  i 

elements  of  Hi,  the  equation  (4.8)  implies  the  Euler  differential  equation  in  D2  (8] 

V FvoiO,  VO,  f,  k)-F0(Q , VO,  z,  k) 

(4.9)  > 

= f A (y,  k )K  (y,  f,  k )E\KO( y,  *)-/?(%*))  ^y, 

■>c2 

and  the  natural  boundary  condition 

(4.10)  Fvo(0,VQ,z-,k)  Uo-  = 0. 

Here,  in  the  arguments  of  Fo  and  F^o.  the  conjugate  variables  O*  and  VO*  are 
omitted  for  simplicity,  and  the  equations  corresponding  to  50*  are  simply  the 
complex  conjugate  of  (4.9)  and  (4.10)  and,  hence,  will  not  be  needed. 

If  no  side  conditions  are  imposed,  the  parametric  function  A should  be  suitably 
chosen  so  that  the  boundary-value  problem  (4.9)  and  (4.10)  has  a unique,  stable 
solution  Oa  in  Hi.  Otherwise,  e.g.  regarding  (4.2)  as  a side  condition,  we  shall 
consider  A to  be  a Langrangian  multiplier.  In  either  case,  Oa  wi"  be  accepted  as  a 
generalized  solution  to  the  ill-posed  problem.  To  be  specific  and  to  demonstrate  the 
flexibility  of  the  proposed  scheme,  we  shall  present  three  examples  which  give 
different  interpretations  of  the  general  principle. 

5.  Examples. 

A.  Continuous  source  with  the  least  noise  intensity.  Perhaps  the  simplest  choice  of 
the  structural  functional  S'  is  the  total  intensity  of  Q defined  as 

(5.1)  ^[Q.O’l^kJ  0(z,k)0*{z,k)dzdk. 

Among  all  solutions  of  (4.2),  we  wish  to  single  out  the  one  that  minimizes  J.  The 
problem  of  minimizing  2T  subject  to  the  constraint  (4.2)  can  be  solved  by  the  method 
of  the  Lagrangian  multiplier  A with  E(u)~  u.  In  this  case  the  Euler’s  equation  (4.9) 
becomes  simply 

(5.2)  a(f,*)=-[  A(y,  k)K(y,  f,  k)dy  = —K*k(z,  k), 

Jo2 

where  K*  designates  the  adjoint  of  K. 

Let  us  introduce  the  integral  operator  It  defined  as 

(5.3)  (y,  k)  = KK*\(y,  k)=  f ( K(y,  z,k)\(y',k)K*(y,  z,k)dz  dy'. 

Jd2  Ju2 

Upon  substituting  (5.2)  into  (4.2),  we  obtain  the  integral  equation  for  A,  in  the 
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operator  notation 

<54)  *A  =-R 

Note  the  £ defined  above  is  a self-adjoint,  Hilbert-Schmidt  operator  on  L2(U 2)  If 
•he  rad.at.on  amplitude  A is  in  the  range  of  * the  Fredholm  integral  equation  ( 4 
has  a unique  solution  given  by  4 

(55>  £—*-'R—(KK*r'R. 

Alternatively  R belongs  to  the  domain  of  the  inverse  ^ l provided  that  [19) 

,fl  ^ I ~t|  [ R(z,  k)dz\  <oo, 

where  W and  {<*„}  are  the  eigenvalues  and  eigenfunctions  of  K respectively 
In  view  of  (5.5),  the  equation  (5.2)  yields  a generalized  (or  pseudo)  solution 

(57)  Ot  = K*(KK*)'R  = TR. 

T l^Tr'V*0  'r  T,o,  the,abL°Ve  reSU,t  T = K*Wr'  with  the  pseudo  inverse 
' ■ ,1  *7  °f  K l9‘'  Whlch  can  be  c°mputed  easily  by  applying  K*  to  the 

equation  (4.2)  and  inverting  the  resulting  equation.  This  pseudo  inverse  has  the 
stamncal  interpretation  as  the  best  linear  estimator  which  minimizes  the  mean-square 

"TS  1 ,Sr.SUrPri,ing  ,hM  ,W°  different  se,ec,ion  ^"Ciples  lead  to  two 
distinct  generalized  solutions  for  a given  image  function  R 

ac  ,hB  C°”‘‘nUOUS  source  with  a Potential.  Next  we  consider  the  source  correlation  Q 

that  theTtr  rUnl  f‘ate  °f  a IiC,,t,OUS  physical  Variable  in  the  continuous  body  D2  such 
that  the  structural  functional  corresponds  to  the  total  “potential  energy” 


(5.8) 


^10,  O*]  - i | {-a(z,  k)VQ  ■ VO*  + b(z,  k )QQ*}  dz  dk. 


where  for  each  k a and I * are  real,  positive,  continuous  functions,  and  a is  continu- 

conclitionTrT??  o ° eqU,'lbrium' ,he  Potential  energy  is  minimum  under  the 
condition  (4.2)  Similar  to  the  previous  example,  one  is  led  to  the  regularized  bound- 
ary value  problem  for  the  Euler  equation 


(5.9) 


(5.10) 


V(aVO)+hO  = -[  A (y,k)K(y,£,k)dy, 
Ju2 


VO  n 


dO 


= 0. 


ho2  r)n\jD, 

den°te  ,he  Neumann  fur|ction  [15J  for  the  above  boundary-value 
problem.  I hen  this  problem  is  equivalent  to 

(5.11)  Od.k).-NK'A(l.k).- [_£*(,, 

where  A is  as  yet  to  be  determined.  To  this  end,  we  use  (5.1 1)  in  (4  2)  to  get  the 
operator  equation  K 

(5I2)  l?*(y,k)~(KNK*)A(y,k)=  R(y,k). 

Again  ^defined  above  is  a self-adjoint,  Hilbert-Schmidt  operator  from  l.2(U2)  into 
itself  If  R satisfies  the  solvability  condition  (5.6),  by  inverting  (5.12),  the  equation 
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(5.11)  gives  the  generalized  solution 

(5.13)  Ch(z,k)  = NK*(KNK*)  ]R(z,k). 

Of  course,  when  a = 0,  b s 2,  we  have  N = /,  an  identity  operator,  and  the  result  (5.7) 
follows  as  a special  case. 

We  remark  that  the  equations  (5.3)  and  (5.12)  for  A are  unstable.  In  actual 
computation,  an  additional  smoothing  procedure,  such  as  Tihonov’s  is  needed  to  ensure 
stability.  The  next  example  will  illustrate  this  procedure. 

C.  Axisymmelric  line  source.  As  shown  in  [20],  for  an  axisymmetric  line  source 
distributed  continuously  along  the  x-axis  in  the  interval  / = [0,  1],  the  integral  equa- 
tion (4.2)  reduces  to 


(5.14) 


KQ(y,k)  = cxp[-ik(ax  - Py)\Q(x,  y,  k)  dx  dy 
7/ 2 

= R(y,k), 


where  z = (x,  y ),  y = (a,  0),  and  a,  /3  are  the  direction  cosines  of  the  lines  of  sight  with 
respect  to  the  x-axis  as  indicated  in  § 2.  Assuming  the  condition  (5.6)  for  the  data  R 
holds,  the  integral  equation  (5.14)  has  a unique  solution,  but  the  solution  does  not 
depend  continuously  on  data.  For  computational  purposes,  the  general  procedure 
outlined  before  can  be  employed  as  a smoothing  scheme.  In  particular,  if  we  choose 
the  residual  function  E(u)=2uu*  and  A = 1/e,  where  e >0  is  a small  parameter,  the 
Euler’s  equation  (4.9)  reads 

(5.15)  e{V  FV0-F0}  = (K*K)Q-K*R, 


and 


(5.16) 


dn  I a/* 


For  instance,  let  us  set  F[Q,  (?*]  = 2{-aV<?  • VO*  + bQQ*}  as  in  the  previous  exam- 
ple. The  equation  (5.15)  takes  the  form 

(5.17)  e{V  (aVQ)  + bQ}  = K*R-(K*K)Q. 

Thereby  we  have  arrived  at  Tihonov’s  regularization  as  a special  case,  where  3 is 
termed  as  the  smoothing  functional.  The  boundary-value  problem  (5.16)  and  (5.17)  is 
well-posed  because  it  is  equivalent  to  the  Fredholm  integral  equation  of  second  kind: 


(5.18)  ( NK*K)Q-eQ  = (NK*)R . 

In  practice  the  smoothed  problem  may  be  solved  by  finite  deferences  or  finite  ele- 
ments methods,  where  the  value  of  e is  usually  chosen  by  trial  and  error.  Numerical 
examples  can  be  found  in  the  paper  [20],  However  a systematic  way  of  choosing  the 
parameter  c was  discussed  in  detail  by  Franklin  [21]  in  one  of  his  interesting  papers  on 
ill-posed  problems. 

Appendix. 

A.l.  Regularization  of  a general  inverse  radiation  problem.  The  regularization 
procedure  presented  in  $ 4 does  not  limit  its  application  to  the  wave  equation.  In  fact 
the  method  applies  to  other  linear  evolution  equations  provided  that  the  far-field 
radiation  data  is  stationary  and  can  be  properly  defined.  Here  we  shall  treat  a linear 
system  of  hyperbolic  equations  of  first  order.  Let  u - (wi,  10.  * ■ • , u„)  be  a vector  field 


A 
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satisfying  the  following  system  of  equations 

( )u  , 

(A.l)  Mu  = — + A(Vx)u  =q(x,t,to),  jreR, 

at 

where  A(x)  is  a matrix-valued,  linear  function  of  .f,  and  q is  a vector-valued,  sta- 
tionary random  function,  describing  the  source  distribution  in  a region  D.  Suppose 
that  (?  denote  the  Riemann  radiation  matrix  [ 1 5]  satisfying  the  adjoint  equation 

(A. 2)  M*d(x-l  t-s)=8(t-s)S(x-i)I, 

where  / signifies  the  unit  matrix.  The  radiation  field  of  (A. 2)  is  then  given  by 

(A. 3)  u(x,  t,  u>)  = Gq(x,  t.  to)  = | | (j(x t - s)q(i,  s,  to)  ds  d£. 

If  < q ) - 0,  we  form  the  covariance  matrices  (or  tensors)  of  u and  q defined  as 
(A. 4)  P(x,  y,  t)  = (u(x,  t + r)u*(y,  r)>,  0(x,  y,  r)=  (q(x,  t + r)q*(y, /)). 

In  view  of  (A. 3)  and  (A. 4),  we  obtain 

(A. 5)  P(x,  y,  t)=  ||  J d(x-lt  + T-s)0(l  fj,  s -s')(j*(y  -i j,  t - s) ds ds'  d£ drj. 

The  above  integral  can  be  shown  to  be  independent  of  t by  the  Fourier  transform  (2.7) 
in  r (setting  c = 1 ) and  by  a change  of  variable 

(A. 6)  P(x,y,k)=[  \G®G](x-iy-rj,k)Q(iv,T)d{df) 

where  (G  0 G)Q  = GQG*  denotes  a tensor  product.  As  in  § 2,  given  the  radiation 
data  on  the  spherical  surface  Sp,  we  define  the  normalized  radiation  covariance  similar 
to  (2.14) 

(A. 7)  R(a,  f3,  k)=  lim  \G  ® <j]  '(t,  y,  k)P(x,  y,  k). 

x.veSp 

P-* 

Corresponding  to  the  kernel  K given  by  (4.1),  we  have 

(A. 8)  K(y,£,Ic)=  lim  [G  ® G)  \x,  y,  k)lG  ® G](x -f,  y - fj,  k ). 

i.yeS„ 

a-* 

Again  we  arrive  at  the  familiar  equation 

(A. 9)  [ K(y,(.k)QU.k)d(  = R(y,k). 

J o 1 

Thus  it  becomes  clear  that  our  regularization  procedure  introduced  in  $ 4 works  in  the 
general  case. 

A.2.  A sl<  thaslic  well-posed  extension.  Without  being  specific,  consider  the 
stochastic  linear  equation 

(A.  10)  /»  = Aq. 

where  A is  bounded  linear  operator  taking  a random  vector  qeHt  into  another 
random  vector  p e //,  where  //.  //,  are  two  complex  Hilbert  spaces.  In  particular  the 
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equation  (A. 10)  may  represent  the  far-field  relation  of  (2.2)  in  the  spectral  (k) 
domain.  Suppose  that  the  random  data  p is  distorted  by  a noise  d in  H.  Instead  of 
(A.  10),  we  must  consider  the  perturbed  equation 

(A.  11)  Aq=p  + d.' 

Assuming  that  all  random  vectors  are  centered,  the  covariances  P,  Q,  D,  R„,  Rpd,  Rqd 
are  defined  in  terms  of  the  inner  products  ( • , -),[•,  • ] in  H and  H i,  respectively. 
For  example 

W,g)=«P,/)(p,g)>,  [Ofu  g.]  = <k  /.Ik.  g.]>, 

(A. 12) 

(Kpqgi. /)  = <(P. /)[<?.  g i ]>■  etc., 
where  /,  g 6 H and  f,,g,  6 Ht. 

We  wish  to  estimate  q by  a bounded  linear  transformation  B of  p,  i.e. 

(A.  13)  q = Bp , 

so  that  the  mean-square  error  is  a minimum 

(A. 14)  e2  = min  <[(q -<?)/>. />]>• 

Following  Franklin  [5],  we  use  (A.  13)  in  (A.  14)  to  determine  B when  the  minimum 
value  e2  is  attained.  By  “completing  the  square”,  it  can  be  shown  that 

(A.  15)  B = RmP'. 

Now  we  make  use  of  (A.  1 1 ) to  get 

(A.  16)  (AOfu  g)  = (Rpqfl,  g)+(/W>.  g). 

or 

(A.  17)  AQ  = RM+Rd„ 

An  elimination  of  Rpd  from  (A.  15)  and  (A.  17)  gives 
(A. 18)  B = (AO  — Rdq)P  ' ■ 

On  the  other  hand  (A.  13)  implies  that 

(A.  19)  0 = BPB*. 

Upon  using  (A.  18)  in  (A.  19)  with  O replaced  by  0,  we  obtain  a nonlinear  well-posed 
equation  for  0: 

(A. 20)  0=(Ad-Rdq)P  '(A0-RJq)\ 

where  Rdq  must  be  suitably  chosen  a priori.  In  contrast  with  Franklin's  approach, 
is  eliminated  from  (A.  15)  to  yield  a linear  equation,  using  (A.  1 1 ) in  terms  of  P and  Rpd. 
which  are  to  be  chosen  judiciously. 
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